You can download these files by clicking on the links
Working with Luminosity
Note that these are large files: Every zip file will range from 355 MB to 411.1 MB
The files are cloud-free composites
In cases where two satellites were collecting data - two composites were produced.
Working with Luminosity
Working with Luminosity
Note that these are large files: Every zip file will range from 355 MB to 411.1 MB
The files are cloud-free composites
In cases where two satellites were collecting data - two composites were produced.
The products are 30 arc second grids, spanning -180 to 180 degrees longitude and -65 to 75 degrees latitude.
Each composite set is named with the satellite and the year (F121995 is from DMSP satellite number F12 for the year 1995)
Working with Luminosity
Three types of images are available:
F1?YYYY_v4b_cf_cvg.tif: Cloud-free coverages tally the total number of observations that went into each 30 arc second grid cell
F1?YYYY_v4b_avg_vis.tif Raw avg_vis contains the average of the visible band digital number values with no further filtering.
F1?YYYY_v4b_stable_lights.avg_vis.tif The cleaned up avg_vis contains the lights from cities, towns, and other sites with persistent lighting, including gas flares.
Working with Luminosity
The DMSP/OLS resolution is approximately 1000m
It comprises six different DMSP satellites:
F10 (1992–1994)
F12 (1994–1999)
F14 (1997–2003)
F15 (2000–2007)
F16 (2004–2009)
F18 (2010–2013)
Working with Luminosity
After downloading, your folder should look like below:
Working with Luminosity
You download all the relevant files for this exercise below:
These are the steps to have the data ready for plotting in ggplot
#Step1: Read country bordersborders <-st_read(dsn="./data/boundaries/all-country-boundaries.shp", quiet =TRUE)#Step2: Creating a folder calledtarget_dir <-"./data/luminosity/F101992"#Step3: Create the target directory if it doesn't existif (!dir.exists(target_dir)) {dir.create(target_dir)}#Step4: Unzip the fileuntar("./data/luminosity/F101992.v4.tar", exdir = target_dir)#Step5: Unzip the file furtherR.utils::gunzip("./data/luminosity/F101992/F101992.v4b_web.avg_vis.tif.gz", remove =FALSE, overwrite=TRUE)#Step6: Read the rasterf1992<-read_stars("./data/luminosity/F101992/F101992.v4b_web.avg_vis.tif")#Step7: Make sure the two files have the same CRSst_crs(f1992)<-st_crs(borders)#Step8: Reduce the resolution of the raster: Higher dx/day values mean lower resolutionf1992b<-st_warp(f1992, st_as_stars(st_bbox(f1992), dx =0.7, dy =0.7))#Step9: Rename the rasternames(f1992b)<-"lights_1992"fig<-ggplot()+geom_stars(data=f1992b)+scale_fill_gradientn(colours=c("black","white"))+#scale_fill_viridis_c(option = "magma",begin = 0.1)+geom_sf(data=borders, linewidth =0.1, fill =NA, color ="white", alpha=0.5)+theme_bw()+theme(panel.grid.major =element_blank(), panel.grid.minor =element_blank())fig
Luminosity in 1992
Luminosity in 1992 in Europe
Raster Values and Properties
The print method gives a summary of the raster’s properties.
print(f1992b)
stars object with 2 dimensions and 1 attribute
attribute(s):
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
lights_1992 0 2 3 3.384738 4 63 515
dimension(s):
from to offset delta refsys x/y
x 1 515 -180 0.7 WGS 84 [x]
y 1 201 75 -0.7 WGS 84 [y]
The class function allows us to see the type of object that we’re dealing with
class(f1992b)
[1] "stars"
Raster Values and Properties
Fundamentally, a stars object is a collection of matrix or array objects.
This is along with additional properties of the dimensions, such as dimension names, Coordinate Reference Systems (CRS).
We can display the structure of a stars object using str.
str(f1992b)
List of 1
$ lights_1992: num [1:515, 1:201] 0 0 0 0 0 0 0 0 0 0 ...
- attr(*, "dimensions")=List of 2
..$ x:List of 7
.. ..$ from : num 1
.. ..$ to : num 515
.. ..$ offset: num -180
.. ..$ delta : num 0.7
.. ..$ refsys:List of 2
.. .. ..$ input: chr "WGS 84"
.. .. ..$ wkt : chr "GEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223"| __truncated__
.. .. ..- attr(*, "class")= chr "crs"
.. ..$ point : logi NA
.. ..$ values: NULL
.. ..- attr(*, "class")= chr "dimension"
..$ y:List of 7
.. ..$ from : num 1
.. ..$ to : num 201
.. ..$ offset: num 75
.. ..$ delta : num -0.7
.. ..$ refsys:List of 2
.. .. ..$ input: chr "WGS 84"
.. .. ..$ wkt : chr "GEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223"| __truncated__
.. .. ..- attr(*, "class")= chr "crs"
.. ..$ point : logi NA
.. ..$ values: NULL
.. ..- attr(*, "class")= chr "dimension"
..- attr(*, "raster")=List of 4
.. ..$ affine : num [1:2] 0 0
.. ..$ dimensions : chr [1:2] "x" "y"
.. ..$ curvilinear: logi FALSE
.. ..$ blocksizes : NULL
.. ..- attr(*, "class")= chr "stars_raster"
..- attr(*, "class")= chr "dimensions"
- attr(*, "class")= chr "stars"
Raster Values and Properties
From this, we see that we have the f1992b object of length 1.
The 2nd line in the str printout specifies the name and contents of the first (and only) item in the list, namely its type, dimensions, and a sample of the first few values.
In our case, the first item is lights_1992, and it is a matrix with 515 rows and 201 columns
Raster Attributes and Values
To access some of the attributes of our objects, we can type:
names(f1992b)
[1] "lights_1992"
We can change names easily in the following wat:
names(f1992b)<-"lum_1992"names(f1992b)
[1] "lum_1992"
Accessing Attributes by Name
We can access our raster’s attributes by using the list access methods
A histogram will help us identify the raster values distribution
#Turning the object into a dataframedf<-st_as_sf(f1992b[1], as_points =FALSE, merge =FALSE)#Creating a histogramggplot(df, aes(x = lum_1992)) +geom_histogram(binwidth =10, color ="white") +ggtitle("Histogram of lum_1992") +xlab("lum_1992") +ylab("Frequency")+theme_bw()
Accessing Raster Values
We thus see that most of values around the world are 0.
This is because 71 percent of the Earth’s surface is water-covered.
The important part is that by transforming the stars object into a dataframe (grid), we can perform the usual operations dedicated to dataframes.
#Step1: Reading the central polygonrelev_poly<-st_read(dsn="./data/rome_center/central_polygon.shp", quiet = T)#Step2: Reading roadsgis_osm_roads <-st_read(dsn="./data/selection.gdb", layer="gis_osm_roads", quiet = T)#Step3: Extracting the bix from the central polygonext2<-st_bbox(relev_poly)#Step4: Ensuring the same crs for the polygon as for the rasterst_crs(relev_poly)<-st_crs(f1992)#Step5: Cropping the rasterrc <-st_crop(f1992, relev_poly)#Step6: Mapping everythingggplot()+geom_stars(data = rc)+geom_sf(data = relev_poly, color="white", fill=NA)+geom_sf(data = gis_osm_roads, color="white", linewidth=0.05)+theme_bw()+scale_fill_gradientn(colours=c("black","white"))
We can easily replace the values inside the raster (grid) by using simple dataframe operations.
For example, let’s replace all the values that have 63 with NA:
Let us now make the NA grids transparent
Writing raster to file
We can of course write the grid as a shapefile for we can turn it back to a raster.
rc3_raster =st_rasterize(rc3)class(rc3_raster)
[1] "stars"
class(rc3)
[1] "sf" "data.frame"
#Step1: Creating a folder calledtarget_dir <-"./data/output"#Step2: Create the target directory if it doesn't existif (!dir.exists(target_dir)) {dir.create(target_dir)}#Step3: Writing the rasterwrite_stars(rc3_raster, './data/output/rc3_raster.tif', update =TRUE)#Step4: Writing the shape filesst_write(rc3, "./data/output/rc3_grid.shp", delete_dsn =TRUE)
Deleting source `./data/output/rc3_grid.shp' using driver `ESRI Shapefile'
Writing layer `rc3_grid' to data source
`./data/output/rc3_grid.shp' using driver `ESRI Shapefile'
Writing 535 features with 1 fields and geometry type Polygon.
Working with More Complex Rasters
Working with lights was simple
In many other contexts, we can have more complex rasters
Temperature data from CEDA is an example of a raster that contains time series
For a start we will simply plot a date to get a better sense of the structure of this type of raster.
library(terra)#Step1: Read the rasterb <-rast("./data/temp/cru_ts4.07.1901.2022.tmp.dat.nc")#Step2: Extract time informationtime_values <-time(b)#Step3: Turning strings into timedate_objects <-as.Date(time_values)#Step4: Select only the dates from the year 1901 and 2022dates_1901 <-unique(date_objects[format(date_objects, "%Y") =="1901"])dates_2022 <-unique(date_objects[format(date_objects, "%Y") =="2022"])#Step5: Find the index of the selected dates in the time_valuesselected_indices_1901 <-which(time_values %in%as.Date(dates_1901))selected_indices_2022 <-which(time_values %in%as.Date(dates_2022))#Step6: Extract raster values for the selected datesselected_rasters_1901 <- b[[selected_indices_1901]]selected_rasters_2022 <- b[[selected_indices_2022]]#Step7: Extract layer nameslayer_names_1901 <-names(selected_rasters_1901)layer_names_2022 <-names(selected_rasters_2022)
Temperature
For a start we will simply plot a date to get a better sense of the structure of this type of raster.
#Step8: Find indices of layers that start with "tmp"#tmp stands for temperature. There are other rasters within these files.indices_tmp_1901 <-grep("^tmp", layer_names_1901)tmp_layers_1901 <- selected_rasters_1901[[indices_tmp_1901]]indices_tmp_2022 <-grep("^tmp", layer_names_2022)tmp_layers_2022 <- selected_rasters_2022[[indices_tmp_2022]]#Step9: Changing raster names to temperaturenames(tmp_layers_1901)<-dates_1901names(tmp_layers_2022)<-dates_2022#Step10: Turning the spatraster into rasterstars_list_1901 <-lapply(tmp_layers_1901, function(raster) st_as_stars(raster::raster(raster)))stars_list_2022 <-lapply(tmp_layers_2022, function(raster) st_as_stars(raster::raster(raster)))#Step11: Use lapply to apply st_apply to each stars object to calculate the meanmean_list_1901 <-lapply(stars_list_1901, function(stars_obj) st_apply(stars_obj, c("x", "y"), mean))mean_list_2022 <-lapply(stars_list_2022, function(stars_obj) st_apply(stars_obj, c("x", "y"), mean))#Step12: Combine the list of results into a single stars objectmean_stars_1901 <-do.call(stars::st_as_stars, mean_list_1901)mean_stars_2022 <-do.call(stars::st_as_stars, mean_list_2022)#Step13: Changing the name of the rastersnames(mean_stars_1901)<-"avg_tmp_1901"names(mean_stars_2022)<-"avg_tmp_2022"
Temperature
Temperature in 1901
fig<-ggplot()+geom_stars(data=mean_stars_1901)+scale_fill_gradient2(low ='blue',high ='red',mid ='yellow', midpoint =5,na.value =NA)+geom_sf(data=borders, linewidth =0.01, fill =NA, color ="white", alpha=0.5)+theme_bw()+theme(panel.grid.major =element_blank(), panel.grid.minor =element_blank())
Now what if we wanted to extract the average temperature for every single year and make a time plot?
Working with Temperature
This is where our programming skills come in handy from earlier lectures
library(terra)# Step 1: Read the rasterb <-rast("./data/temp/cru_ts4.07.1901.2022.tmp.dat.nc")# Step 2: Extract time informationtime_values <-time(b)# Step 3: Turning strings into timedate_objects <-as.Date(time_values)# Step 4: Get unique yearsunique_years <-unique(format(date_objects, "%Y"))#All years take too long, so we will only choose two yearsunique_years <-c("2021", "2022")# Step 5: Create an empty list to store mean stars for each yearmean_stars_list <-list()
Working with Temperature
This is where our programming skills come in handy from earlier lectures
# Step 6: Iterate over each unique yearfor (year in unique_years) {cat("Calculating average for year", year, "\n")# Select dates for the current year selected_dates <-unique(date_objects[format(date_objects, "%Y") == year])# Find the index of the selected dates in the time_values selected_indices <-which(time_values %in%as.Date(selected_dates))# Extract raster values for the selected dates selected_rasters <- b[[selected_indices]]# Extract layer names layer_names <-names(selected_rasters)# Find indices of layers that start with "tmp" indices_tmp <-grep("^tmp", layer_names)# Extract temperature layers tmp_layers <- selected_rasters[[indices_tmp]]# Change raster names to temperaturenames(tmp_layers) <- selected_dates# Convert spatraster into raster stars_list <-lapply(tmp_layers, function(raster) st_as_stars(raster::raster(raster)))# Use lapply to apply st_apply to each stars object to calculate the mean mean_list <-lapply(stars_list, function(stars_obj) st_apply(stars_obj, c("x", "y"), mean))# Combine the list of results into a single stars object mean_stars <-do.call(stars::st_as_stars, mean_list)# Change the name of the rasternames(mean_stars) <-paste("avg_tmp", year, sep ="_")# Add the mean stars to the list mean_stars_list[[year]] <- mean_starscat("Average calculation for year", year, "completed\n\n")}
Calculating average for year 2021
Average calculation for year 2021 completed
Calculating average for year 2022
Average calculation for year 2022 completed
Working with Temperature
This is where our programming skills come in handy from earlier lectures
# Combine all the mean stars into a single stars object# Initialize an empty data frameresult_df <-data.frame()# Loop through each year and calculate the meanfor (year innames(mean_stars_list)) {#Creating an object where we look through the years avg_tmp_col_name <-paste0("avg_tmp_", year)#Calculating mean and keeping track of the year mean_value <-mean(mean_stars_list[[year]][[avg_tmp_col_name]], na.rm =TRUE)#Saving everything into a dataframe result_df <-rbind(result_df, data.frame(year =as.numeric(year), average_temperature = mean_value))}# Writing the df to a dataframewrite.csv(result_df, "./data/output/average_temp_2021_2022.csv", row.names =FALSE)
Working with Temperature
So, we wrote our dataframe to a csv file.
Note that we did this dataframe for two years only: 2021 and 2022.
To complete the loop for all years, put a # in front of unique_years <-c("2021", "2022")
This means that unique_years is an object made out of unique(format(date_objects, "%Y"))
In other words, all years from 1901 until 2022. But this takes a long time (i.e. a few hours): multiply how long it took for 2021 and 2022 by 50.
To save you time, I have already done all the years on my machine. You can download the csv file with all the years.
Put it in the relevant folder
Working with Temperature
#write.csv(result_df, "./data/output/average_temp_1901_2022.csv", row.names = FALSE)#Reading the fileaverage_temp<-read.csv("./data/output/average_temp_1901_2022.csv")#Subsetting to below 2001average_temp_1901_2000<-subset(average_temp, year<2001)#Calculating the 1901-2000 averageaverage_temp_1901_2000<-mean(average_temp_1901_2000$average_temperature)#Calculation the deviations from the 1901-2000 averageaverage_temp$dev_1901_2000_avg<-average_temp$average_temperature-average_temp_1901_2000
Working with Temperature
Creating a time trend
#Mapping the Temperature Deviaitonsover timeggplot(average_temp, aes(x = year, y = dev_1901_2000_avg)) +geom_line() +geom_point() +theme_bw()+labs(title ="Time Trend of Average Temperature",x ="Year",y ="Average Temperature")
Working with Temperature
Creating a time trend
Working with Temperature
Creating a bar plot with the values
# Create a new column for bar color based on the conditionaverage_temp$bar_color <-ifelse(average_temp$dev_1901_2000_avg >0, "red", "blue")# Create the bar plotggplot(average_temp, aes(x = year, y = dev_1901_2000_avg, fill = bar_color)) +geom_bar(stat ="identity", color ="black") +scale_fill_manual(values =c("blue", "red"), guide="none") +theme_bw() +labs(title ="Time Trend of Average Temperature",x ="Year",y ="Difference from 1901-2000 avearge (degrees C)")
Working with Temperature
Creating a bar plot with the values
Working with Temperature
This is very close to the official figures produced by NOAA (National Ceneter for Environmental Information)
Without realizing, we have already covered an important component of working with rasters: raster algebra.
Raster algebra is the application of functions on a cell-by-cell basis, repeatedly for all pixels of one or more overlapping rasters.
Example:
Raster Algebra
Without realizing, we have already covered an important component of working with rasters: raster algebra.
Raster algebra is the application of functions on a cell-by-cell basis, repeatedly for all pixels of one or more overlapping rasters.
Example:
Raster Algebra
Without realizing, we have already covered an important component of working with rasters: raster algebra.
Raster algebra is the application of functions on a cell-by-cell basis, repeatedly for all pixels of one or more overlapping rasters.
Example:
Raster Algebra
Without realizing, we have already covered an important component of working with rasters: raster algebra.
Raster algebra is the application of functions on a cell-by-cell basis, repeatedly for all pixels of one or more overlapping rasters.
Example:
Raster Algebra
Without realizing, we have already covered an important component of working with rasters: raster algebra.
Raster algebra is the application of functions on a cell-by-cell basis, repeatedly for all pixels of one or more overlapping rasters.
Example:
Raster Algebra
Without realizing, we have already covered an important component of working with rasters: raster algebra.
Raster algebra is the application of functions on a cell-by-cell basis, repeatedly for all pixels of one or more overlapping rasters.
Example:
Raster Algebra
Raster Algebra Operations
Arithmetic: +, -, *, /
Logical: <, <=, >, >=, ==, !=, !
We can also apply the following operations to individual rasters: is.na, abs, round, sqrt, log, etc.
Raster Algebra Illustration
Let’s focus once again on Italy
#Raster Addition
We can add the two rasters together in the following way:
#Raster Average
But it probably makes sense to create an average of the two rasters
#Raster Average
As we showed previously, can make raster raplacements by using the dataframe operations
More Raster Algebra
As we saw previously, satellite luminosity is used as a proxy for economic activity
Luminosity in 2013
Working with Luminosity
These are the steps to have the data ready for plotting in ggplot
#Step1: Read country bordersborders <-st_read(dsn="./data/boundaries/all-country-boundaries.shp", quiet =TRUE)#Step2: Creating a folder calledtarget_dir <-"./data/luminosity/F101992"#Step3: Create the target directory if it doesn't existif (!dir.exists(target_dir)) {dir.create(target_dir)}#Step4: Unzip the fileuntar("./data/luminosity/F101992.v4.tar", exdir = target_dir)#Step5: Unzip the file furtherR.utils::gunzip("./data/luminosity/F101992/F101992.v4b_web.avg_vis.tif.gz", remove =FALSE, overwrite=TRUE)#Step6: Read the rasterf1992<-read_stars("./data/luminosity/F101992/F101992.v4b_web.avg_vis.tif")#Step7: Make sure the two files have the same CRSst_crs(f1992)<-st_crs(borders)#Step8: Reduce the resolution of the rasterf1992b<-st_warp(f1992, st_as_stars(st_bbox(f1992), dx =0.7))#Step9: Rename the rasternames(f1992b)<-"lights_1992"fig<-ggplot()+geom_stars(data=f1992b)+scale_fill_gradientn(colours=c("black","white"))+#scale_fill_viridis_c(option = "magma",begin = 0.1)+geom_sf(data=borders, linewidth =0.1, fill =NA, color ="white", alpha=0.5)+theme_bw()+theme(panel.grid.major =element_blank(), panel.grid.minor =element_blank())fig
Luminosity in 1992
Luminosity in 1992 in Europe
library(purrr)#Identifying coordinates of Northern pointsnorway<-subset(borders, SOVEREIGNT =="Norway")nc_geom <-st_geometry(norway)list_points_norway<-nc_geom[[1]]list_points_norway<-flatten(list_points_norway)list_points_norway2<-as.data.frame(do.call(rbind, list_points_norway))names(list_points_norway2)<-c("lon_x", "lat_y")max_lat_y_eu<-max(list_points_norway2$lat_y)#Identifying coordinates of Southern pointsgreece<-subset(borders, SOVEREIGNT =="Greece")nc_geom <-st_geometry(greece)list_points_greece<-nc_geom[[1]]list_points_greece<-flatten(list_points_greece)list_points_greece2<-as.data.frame(do.call(rbind, list_points_greece))names(list_points_greece2)<-c("lon_x", "lat_y")min_lat_y_eu<-min(list_points_greece2$lat_y)#Identifying coordinates of Eastern pointsukraine<-subset(borders, SOVEREIGNT =="Ukraine")nc_geom <-st_geometry(ukraine)list_points_ukraine<-nc_geom[[1]]list_points_ukraine<-flatten(list_points_ukraine)list_points_ukraine2<-as.data.frame(do.call(rbind, list_points_ukraine))names(list_points_ukraine2)<-c("lon_x", "lat_y")max_lon_x_eu<-max(list_points_ukraine2$lon_x)#Identifying coordinates of Western pointsportugal<-subset(borders, SOVEREIGNT =="Portugal")nc_geom <-st_geometry(portugal)list_points_portugal<-nc_geom[[1]]list_points_portugal<-flatten(list_points_portugal)list_points_portugal2<-as.data.frame(do.call(rbind, list_points_portugal))names(list_points_portugal2)<-c("lon_x", "lat_y")min_lon_x_eu<-min(list_points_portugal2$lon_x)
Luminosity in 1992 in Europe
Mapping
fig<-ggplot()+geom_stars(data=f1992b)+scale_fill_gradientn(colours=c("black","white"))+geom_sf(data=borders, linewidth =0.1, fill =NA, color ="white", alpha=0.5)+coord_sf(xlim =c(min_lon_x_eu-3, max_lon_x_eu+3), ylim =c(min_lat_y_eu-1, max_lat_y_eu-6), expand =FALSE)+theme_bw()+theme(panel.grid.major =element_blank(), panel.grid.minor =element_blank())fig
Luminosity in 1992 in Europe
Mapping
Luminosity in 1992 in the US
We can do the same about the US
min_lon_x_us<-(-130)max_lon_x_us<-(-57)min_lat_y_us<-(26)max_lat_y_us<-(53)fig<-ggplot()+geom_stars(data=f1992b)+scale_fill_gradientn(colours=c("black","white"))+geom_sf(data=borders, linewidth =0.1, fill =NA, color ="white", alpha=0.5)+coord_sf(xlim =c(min_lon_x_us-3, max_lon_x_us+3), ylim =c(min_lat_y_us-3, max_lat_y_us+3), expand =FALSE)+theme_bw()+theme(panel.grid.major =element_blank(), panel.grid.minor =element_blank())fig
Luminosity in 1992 in the US
We can do the same about the US
Luminosity in 2013
Luminosity in 2013 in Europe
Luminosity in 2013 in the US
Comparing Luminosity
Comparing Luminosity
Positive Changes: 2013-1992>0
Comparing Luminosity
Here is how we do that:
#Step8: Reduce the resolution of the raster: Higher dx/day values mean lower resolutionf1992c<-st_warp(f1992, st_as_stars(st_bbox(f1992), dx =0.5, dy =0.5))#Step9: Rename the rasternames(f1992c)<-"lights_1992"#Step8: Reduce the resolution of the rasterf2013c<-st_warp(f2013, st_as_stars(st_bbox(f2013), dx =0.5, dy =0.5))#Step9: Rename the rasternames(f2013c)<-"lights_2013"#Step10: Turning to gridf1992bpo<-st_as_sf(f1992c, as_points =FALSE, merge =FALSE)f2013bpo<-st_as_sf(f2013c, as_points =FALSE, merge =FALSE)#Step11: Spatially intersection 1992 and 2013joined_data <-st_join(f1992bpo, f2013bpo, join = st_intersects)#Switched off spherical geometry (s2)sf::sf_use_s2(FALSE)#Step12: Keeping only the ADMIN variableborders2<-subset(borders, select =c("ADMIN"))#Step13: Simplifying the world polygonsborders2<-st_simplify(borders2, dTolerance =0.005)#Step14: Spatially intersecting borders and gridsjoined_data_w<-st_join(joined_data, borders2, join = st_intersects)#Step15: Removing seasjoined_data_w2<-subset(joined_data_w, !is.na(ADMIN))#Step16: Calculating changes in luminosityjoined_data_w2$diff_lum<-joined_data_w2$lights_2013-joined_data_w2$lights_1992joined_data2<-subset(joined_data_w2, diff_lum>0)joined_data3<-subset(joined_data_w2, diff_lum<(-5))
Comparing Luminosity
Here is how we do that:
Positive Changes: 2013-1992>0
ggplot()+geom_sf(data=joined_data2, color=NA, aes(fill=diff_lum))+scale_fill_gradientn(colours=c("yellow", "red"), na.value =NA)+geom_sf(data=borders, linewidth =0.01, fill =NA, color ="black", alpha=0.5)+theme_bw()+theme(panel.grid.major =element_blank(), panel.grid.minor =element_blank())
library(terra)library(sf)#Step1: Reading the rasterb <-rast("./data/temp/cru_ts4.07.1901.2022.tmp.dat.nc")#Step2: Turning it into stars objectstars_obj1 <-st_as_stars(raster::raster(b$tmp_1424))#Step3: Renaming the rasternames(stars_obj1)<-time(b$tmp_1424)#Step4: Ensuring the the same CRSst_crs(stars_obj1)<-st_crs(border)
Temperature
Temperature in 2019
fig<-ggplot()+geom_stars(data=stars_obj1)+scale_fill_gradient2(low ='blue',high ='red',mid ='yellow', midpoint =5,na.value =NA)+geom_sf(data=borders, linewidth =0.01, fill =NA, color ="white", alpha=0.5)+theme_bw()+theme(panel.grid.major =element_blank(), panel.grid.minor =element_blank())
Temperature
Temperature in 2019
Temperature Difference
#Step1: Reduce the resolution of the raster: Higher dx/day values mean lower resolutionstars_1901<-st_warp(mean_stars_1901, st_as_stars(st_bbox(mean_stars_1901), dx =1, dy =1))stars_2022<-st_warp(mean_stars_2022, st_as_stars(st_bbox(mean_stars_2022), dx =1, dy =1))#Step2: Turning to gridgrid1901<-st_as_sf(stars_1901, as_points =FALSE, merge =FALSE)grid2022<-st_as_sf(stars_2022, as_points =FALSE, merge =FALSE)#Step3: Spatially intersection 1992 and 2013joined_data <-st_join(grid1901, grid2022, join = st_intersects)#Step4: Calculating the differencejoined_data$diff_temp<-joined_data$avg_tmp_2022-joined_data$avg_tmp_1901joined_data2<-subset(joined_data, diff_temp>4)joined_data3<-subset(joined_data, diff_temp<(-5))
Plotting Changes
Temperature difference between 2022 and 1901
fig<-ggplot()+geom_sf(data=joined_data, color=NA, aes(fill=diff_temp))+scale_fill_gradientn(colours=c("yellow", "blue"), na.value =NA)+geom_sf(data=borders, linewidth =0.01, fill =NA, color ="black", alpha=0.5)+theme_bw()+theme(panel.grid.major =element_blank(), panel.grid.minor =element_blank())
Plotting Changes
Temperature difference between 2022 and 1901
Plotting Changes
Temperature difference between 2022 and 1901 > 4
fig<-ggplot()+geom_sf(data=joined_data2, color=NA, aes(fill=diff_temp))+scale_fill_gradientn(colours=c("yellow", "blue"), na.value =NA)+geom_sf(data=borders, linewidth =0.01, fill =NA, color ="black", alpha=0.5)+theme_bw()+theme(panel.grid.major =element_blank(), panel.grid.minor =element_blank())
Plotting Changes
Temperature difference between 2022 and 1901 > 4
Conclusion
We have covered a lot of operations related to rasters